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The microcanonical properties of a two dimensional system of A*' classical particles interacting 
via a smoothed Newtonian potential as a function of the total energy E and the total angular 
momentum L are discussed. The two first moments of the distribution of the linear momentum of a 
given particle at a fixed position show that, (a) in average the system rotates like a solid body (b) the 
equipartition theorem has to be corrected by a term proportional to the square of the fiuctuations of 
the inertial momentum of the system. In order to estimate suitable observables a numerical method 
based on an importance sampling algorithm is presented. The entropy surface S shows a negative 
specific heat region at fixed L for all L. Observables probing the average mass distribution are used 
to understand the link between thermostatistical properties and the spatial distribution of particles. 
In order to define a phase in non-extensive system we introduce a more general observable than the 
one proposed by Gross and Votyakov [Eur. Phys. J. B 15, 115 (2000)]. This observable is the sign 
of the largest eigenvalue of the entropy surface curvature. If it is negative then the system is in a 
pure phase; if it is positive then the system undergoes a first order phase transition. At large E the 
gravitational system is in a homogeneous gas phase. At low E there are several collapse phases; at 
L = there is a single cluster phase and for L 7^ there are several phases with two clusters, the 
relative size of the clusters depends on L. All these pure phases are separated by first order phase 
transition regions. The signal of critical behaviour emerges at different points of the parameter space 
{E,L). We also discuss the ensemble introduced in a recent pre-print by Klinko and Miller; this 
ensemble is the canonical analogue of the one at constant energy and constant angular momentum, 
i.e. it is defined at constant temperature and constant angular velocity. We show that a huge loss 
of informations appears if we treat the system as a function of intensive parameters: besides the 
known non-equivalence at first order phase transitions, there exit in the microcanonical ensemble 
some values of the temperature and the angular velocity for which the corresponding canonical 
ensemble does not exist, i.e. the partition sum diverges. 



I. INTRODUCTION 



The thermostatistical properties of systems of N classical particles under a long-range attractive potential have 
been extensively studied since the seminal work of Antonov One of their more specific and interesting 

properties is that they are unstable for all N ^ and therefore not thermodynamically extensive, i.e. they exhibit 
negative specific heat regions even when the system is composed by a very large number of particles. 

It is quite natural to ask whether the total angular momentum L, which is an integral of motion for systems 
of relevance in the astrophysical context, plays a non-trivial role on the equilibrium properties of these systems. 
Indeed L is considered as an important parameter in order to understand the physics of systems like galaxies [p|-p^ ; 
globular clusters [pT|~p^ ; molecular clouds in multifragmentation regime [ p^jp^ which might eventually lead to stellar 
formation |l8|-|2^]. 

Previous works have already studied the effect of L in the mean field limit with a simplified potential and imposing 
a spherical symmetry or atL = 

Our work, presented in this paper, is an attempt to overcome some of these approximations. 

Thermodynamical equilibrium does not exist for Newtonian self-gravitating systems, due both to evaporation of 
stars (the systems are not self-bounded) and short distance singularity in the interaction potential. However there 
exists intermediate stages where these two effects might be neglected and a quasi-equilibrium state might be reached 
(dynamical issues like ergodicity, mixing or "approach to equilibrium" ||^,|2^,^ are not considered in this paper). In 
order to make the existence of equilibrium configurations possible we have, first, to bound the system in an artificial 
box and, second, to add a short distance cutoff to the potential. The latter point can be seen as an attempt to take 
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into account the appearance of new physics at very short distances (about the influence of this short cut see p8| , p9| ). 
Another way to avoid the difficulties due to the short distance singularity is to describe the function of distribution 
of the "stars" within a Fermi-Dirac statistic ■ 

The box breaks the translational symmetry of the system, therefore the total linear momentum P and angular 
momentum L are not conserved. Nevertheless we assume that the equilibration time is smaller than the characteristic 
time after which the box plays a significant role Therefore P and L are considered as (quasi-)conserved 

quantities. We put the center of the box at the center of mass Rcm which is also set to be the center of the 
coordinates. Therefore P = 0. 

As already mentioned, self-gravitating systems are non-extensive and a statistical description based on their in- 
tensive parameters should be taken with caution since the different statistical ensembles are only equivalent at the 
thermodynamical limit far from first order transitions (see Sec. HID). Moreover, this limit is required in order to 
define phases and phase transitions if one fixes the intensive parameters |31 . In contrast, the microcanonical ensemble 



(ME) does not require this limit, it allows a classification of phase transitions for finite size systems [^,^. Hence 
the considered system is studied within the natural ME framework. 

In order to perform the computation in a reasonable time we have to consider a two dimensional system. 

Th e pap er is organized as follow: In section ^ we recall the analytical expressions for entropy and its derivatives 
(Sec. [I A ), clarify the definition of phase transitions for non-extensive systems proposed in P2| (Sec. II B ), d iscus s 
the two first moments of the distribution of the linear momentum of a given particle at a fixed position (Sec. |II ( 



and present a numerical method based on an importance sampling algorithm in order to estimate suitable observables 
(Sec. [ID). Numerical results are presente d in s ection III ; the link between the average mass distribution and the 
thermostatistical properties is made in Sec. [II B, In Sec. [II C we use the definition of phase introduced in Sec. [IB 



to draw the phase diagram of the self-gravitating system as a function of its energy E and angular momentum L. 
We discuss the ensemble introduced in jl^ ; therein this ensemble is used to treat another model of rotating and 
self-gravitating system. This ensemble is a function of the (intensive) variables conju gate o f E and L. For our model 
we show how the predictions using this ensemble are inaccurate and misleading (Sec. HI D| ) . Results are summarized 



and discussed in section [V 



II. MICROCANONICAL PROPERTIES 



A. Microcanonical definitions 



Consider a system of N classical particles on a disk of radius R whose interaction is described by a Plummer 
softened potential P^,p5[ 



GmiTrij 



(1) 



where and = {qj, qf} are the mass and position of particle i respectively, s is the softening length and G is the 
gravitational constant. The fixed total energy E is described by the Hamiltonian 



PL 

2m, 



(2) 



where Pj = {pl,pf} is the linear momentum of particle i, f = "^Zkj Vij ^^id g is a 2A''-dimensional vector whose 
coordinates are {Qi, ■ ■ ■ , qj^}, representing the spatial configuration and is an element of the spatial configuration 
space Vc, qeVcC M^^. 

The entropy S is given through the Boltzmann's principle (the Boltzmann constant is set to I) 

SiE,L,N)^log{W{E,L,N)), (3) 

where W{E, L, N) is the volume of the accessible phase-space at E, L and N fixed (under the assumptions given in 



(4) 
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where x = qjpf — qfpj . After integration over the momenta Eq. (Q) becomes p5| , ^ 



W{E,L,N)^c[ dq^E^-'^l'^, (5) 

where C = (277/1)^^ m )r(jV-3/2) ^ constant, / = 22i is the inertial momentum and Er = E — ^ — (p 

the remaining energy . From the point of view of the remaining energy, if L 7^ we can already notice that the 
equihbrium properties are the results of a competition between two terms; the rotational energy -^j and the potential 
energy ip. The former tries to drive the particles away from the center of mass in order to increase / whereas the 
latter tries to group the particles together in order to decrease ip, but since the center of mass is fixed this will lead 
to a concentration of particles near the center and consequently will decrease /. 
The microcanonical temperature T is defined by 

1 ^ dS ,N-5/2, 



where (•) is the microcanonical average 



The angular velocity ui is defined as minus the conjugate force of L times T [p7[ 



We also define 7/? has the conjugate of L 



2. 



uj = 2L7. (10) 



B. Phase and phase transitions 



For small or self-gravitating systems a phase (or a phase transition) can not be defined in the usual way using 
for example the Lee and Yang singularities pl| ] since these singularities show up only at the thermodynamical limit. 
Invoking the thermodynamical limit when studying small systems washses out all the finite size effects that may 
lead to new phenomena, (e.g. isomerisation of metallic clusters [Q, multifragmentation of nuclei [^) and for self- 
gravitating systems the thermodynamical limit does not exist. Hereafter, we will call "Small systems", those where 
the range of the forces is of the order of the system size (e.g. metallic clusters, nuclei and self-gravitating systems) 
and also systems without proper thermodynamical limit (e.g. unstable systems pO|). 

In a recent paper [ p2| definitions for pure phases and phase transitions (first and second kind) based on the local 
topology of the microcanonical entropy surface has been proposed. In the following we first fix some notations and 
then recall the definitions. 

Consider the microcanonical ensemble (ME) of an isolated physical system. Its associated entropy S{X) is a 
function of J\f "extensive" dynamical conserved quantities X = {X^, . . . ,X^}. Note that X may not contain all 
the dynamical conserved quantities and for simplicity all these parameters are considered as being continuous. The 
Jacobian ofS{X = Xq) is noted by Js(Xo) = II ox^dxi H^oi its eigenvalues are {Ai, . . . , Aa^} where Ai > A2 > . . . > Xf^ 
and the determinant of Js(X) is Ds = AI • • • Xj^. 

In phase transitions are defined "&?/ the points and regions of non-negative curvature of the entropy surface 
[. . . ] as a function of the mechanical quantities" . Therein the sign of Ds is put forward as a measure the concavity 
of S (its negative curvature) so that at first order phase transition 

sign (I?5)= sign ((-l)-^+i). (11) 

Though this condition is necessary it is not sufficient in the general case. In fact 5 is a non-concave function at Xq if 
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Ai >0, 



(12) 



i.e. if at least one eigenvalue of J5 is non-negative. Note that in the two dimensional example model studied in |3S 



in order to illustrate the definition (11), A2 is always negative therefore the sign of Ds is simply minus the one of Ai, 
and the conditions (|ll|) and ([l^) are equivalent. 

By using sign(Ai), one can somewhat extend or clarify the classification of phase transitions in non-extensive 
systems when the entropy is a function of > 1 variables for any set of the parameters Xq 

• A single pure phase if \i{Xq) < 0. 



A first order phase transition if Xi{Xo) > 0. As mentioned in [ p2y4l| the depth of the associated entropy intruder 
is a measure of the intra-phase surface tension; how to define and to measure these depth when J\f > I will be 
discussed elsewhere. Note that some eigenvalues might still be negative although Xi{Xq) > just like in the 
model presented in |32| . In this case "good" order parameters are linear combinations of the eigenvectors whose 
eigenvalues are positive. 

If Xi{Xo) — and Ai is the only zero eigenvalue and VvjAi — 0, where vi is the eigenvector associated with Ai 
then there is a second order phase transition at Xa- 

If several eigenvalues obey A^ = and Vv; A^ — for i = 1, . . . , n < A/" then Xq is a multicritical point. 



C. Momentum average and dispersion 

In this section we derive the average and the dispersion of the linear momentum of a particle, we also compute its 
mean angular velocity and relate it to the one of the system as defined in Eq. (^. 

The derivation of {Pklq^ the average momentum of particle k at fixed position q^. (while the other particles are 
free) is similar to that of W. Details of the derivation can be found in Appendix and the result is 

2 

^Pk)q^ = L{I^^)q^mk ^ e^sqiea, (13) 

a.S=l 

where e is the antisymmetric tensor of rank 2 and the unit vector of coordinate a. Equation dl^ ) shows that {Pk)q^. 
is a vector perpendicular to q^, whose module is a function of ||qj,||. In other words the orbit of a particle is in the 
mean circular (this result is expected since the system is rotationally symmetric). One can compute (wfc)q^ the mean 
angular velocity of k at distance \\q}^\ by first considering {Lk)q^ the mean angular momentum of k at distance 

{Lk)q, ^qk>^ {Pk)q, 

= L{r')qJk, (14) 

where Ik — rakq^- The angular mean velocity of a particle on a circular orbit is classically linked to {Lk)q^. by 

{Lk)q^ ^ {ujk)qjk- (15) 



We can identify {ujk)q^ in ( p^ ) as 



(^fe)g, =i(Og,- (16) 



The dependence of (w)q^ on Hq^H is of the order 1/N {{I ^)q^ = {I ^) + 0{N ^)), therefore for large N we can 
write (see Eq. (||)) 

{L0k)q^^L{r')^{L0). (17) 

For large N the mean angular velocity is the same for all the particles at any distance from the center, in other 
words the system in the mean rotates like a solid body. Moreover {ijjk)q^ corresponds to the thermostatistical angular 
velocity u defined by Eq. (^). These are also classical results for extensive systems at low L [^,^. Note also that 

these results do not depend on the interaction potential Lp. 

The momentum dispersion cjp^ can also be derived. Using Eq. ( p^ ) and Eq. (|A5|), we get for large N 
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^2^+hL'm,{{I-')-{rr) (18) 



The second term of Eq. (18) is proportional to the square of the dispersion of and to {Ik — mkql). When this 
term vanishes relatively to the first one, e.g. when the fluctuations of I^^ are small, or at high energy (low /3) and 
low L, we recover the usual dispersion of the Maxwell-Boltzmann distribution. This term also gives a correction to 
the usual equipartition theorem; for large N 

at 



- T + {{I-') - (i-r) . (19) 

where (Ek) is the average internal kinetic energy (without the contribution from the collective rotational movement) 
of particle k. Again this correction is position-dependent via Ik- In the regimes where the fluctuations of the mass 
distribution can not be neglected in Eqs. (|l8| ) and (|l^) an estimation of the temperature based on the velocity 
dispersion would show that the temperature is smaller in the core than at the edge. 



D. Numerical method 

From now on we set rui = m, Wi ^ 1, ■ ■ ■ , N and use the following dimensionless variables: (i) E ^ e = q^^^^-i ; 
(ii) L P = = 2gSrn-^ ' (iii) * ^ ^ = (i^) 9 ^ ^ = (v) V'^ ^ Va\ (vi) = m^^f = 

7^ ^i<3 ^^2+i^l^_r^)2 ' (vii) I^I'^Ei'rf- The weight is now 

W{e,n)=C' [ (20) 

Jv^ \I 

where = e — j — is the dimensionles remaining energy and C a constant. Later on we no longer write this 
constant since it plays no significant role (it only shifts the entropy by logC). The derivatives of entropy (/3, w, . . . ) 
are now dimensionless quantities. 

One usually estimates ( ^ ) by means of some Monte Carlo algorithm, updating the positions q by some small 
amount (5q in order to get a good pass acceptance and using the configuration weight W(r) = -^f-r in the 
Metropolis pass. Unfortunately this strategy does not really work (within a reasonable CPU-time), because the 2N- 
dim configuration weight-landscape at fixed e and Q, presents troughs and high peaks p3[ , so exploring the total 
configuration-space (or at least a subset containing the highest peaks) would take a very long, in practice infinite, 
time. This weight-landscape looks like the energy-landscape found in spin-glass systems. 

The strategy we have adopted is described in the following. First we can rewrite (po|) 



W (e, n)^ Jdl dc^ Bg (/, 0) -^j^r~"\ (21) 



where Bg {I, (f) — dr5 {I'{r) — 1)5 {(j)'{r) — (j)); Bg (/, (j>) is the density of spatial configurations at given / and 
Given Bg we can compute W, S and its derivatives for any e and e.g. 

- 1^ 

_ TV - 5/2 / dl dc^ Bg (/, (j)) /-3/2e^- 
The expectation value (O) of an observable 0{r) can be estimated if we know Bg and {O) ^ ^ 



(22) 



Jdldct> {0),^^Bg{IA)I-"'er"' 
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where 



ldrO{r)5{I' -I)5{(l>' 
Bg{IA) 



(24) 



Now, we have to compute Bg{I,(j)) and {O) j ^. A priori Bg{I,(j)) is highly peaked around the values of / and 
that describe the gas (disordered) phase and should drop down to the edges. Nevertheless we need a good estimate 
of Bg(I,(j3) for almost all values taken by (/, 0) even when Bg[I,(f>) is very small comparing to its maximum. For 

example at small total energy e only the part of Bg{I, 0) for which er = e — ^ 0>O will contribute to the 

integral (^). 

In order to get a good estimate of Bg we used an iterative scheme inspired by multicanonical algorithms [^-^. 
The usual multicanonical task is to compute the free energy as a function of the total energy. Here we have to compute 
Bg as a function of / and (p. The updating scheme presented in |^] was used; one of the reason for this choice is that 
although it has been given for a one dimensional task it can be trivially extended to bi-variate problem. Furthermore 
we added a blocking mechanism: once we estimate that enough information has been collected on a given region of the 
parametric space (/, 4>) then we tag it as "locked" so that it will not be visited during next iterations. This mechanism 
enables the program to spread more quickly over the parametric space and save computation time comparing to usual 
multicanonical algorithms. Details on this blocking mechanism will be given elsewhere H|. 



150 



^ 100 

II 



50 



10 



i=l 
i=4 
i=7 
i=10 



(a) 



-4^ 



(b) 



FIG. 1. Estimate of the density of state Bg (a) and histogram of the visited states C (b) for 7 = 3 at different iteration steps 
i of the multicanonical algorithm as a function of the potential energy Panel (a) shows how Bg is built step by step; Bg is 
an extremely peaked function, the log of the ratio between its maximum and its minimum is about 120. Without the blocking 
mechanism (see text) d would have been non null for all values of (j> visited during previous steps j < i- In panel (b) we see 



that the algorithm does no longer visit "well-known" regions 



> 



-1.5) already after four steps. 



In the present paper we present results for N — 20, and cr — 0.05. Figure |l] shows a slice of Bg{I, (j)) for / = 3 at 
different iteration steps i (Fig. [l].a). We have also plotted the histogram C{I — 3,0) of the visited region in order 
to illustrate the blocking mechanism (Fig. |l|.b). As expected Bg is strongly peaked around the disordered region 
(j) « —1 (this value correspond to the mean of (f) over randomly generated spatial configurations); after 10 iterations 
the ratio between the maximum and the minimum of Bg is ~ exp 120. This ratio increases exponentially with N, 
e.g. at = 10 its value is « exp 80. This is the main reason why we could not study systems with larger N within 
the current algorithm, but we should add that no qualitative changes have been noted between = 10 and A^ = 20 
(some preliminary studies support this remark for N = 50). 
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III. RESULTS 



A. Entropy and its derivatives 




FIG. 2. Entropy surface S (e, Q,), the mesh lines are at constant e or constant Q,. The thick line is close to the pojection of 
the T = {S ^ ~oo) isotherm. A convex intruder at constant Q and for a certain energy range (e.g. —2 < e < for Q = 20) 
can be seen for all fi. 5* is not defined in the forbidden region; there the remaining energy is negative for all q G Hc. 

Figure |^ shows the entropy surface 5 as a function of e and 17. The ground state energy eg(il) (thick line in Fig. |^) 



increases with fl; eg classically corresponds to = implying S = 



-oo. 



For all Q, eg{Q) is a concave function of Q, 



i.e. < 0; at high il ( fl > 12) it is almo st line ar — > 0^. These properties show that the set {e, il} over which 
S is defined is not convex Bfl], and in Sec. HID we discuss some consequences resulting from a non-convex domain 
of definition. 

At fixed f2, ^(e) is not concave for all e but shows for some energy interval a convex intruder which signals a first 

order phase transition with negative specific heat < 0^ |33|. This can be better viewed by plotting /3(e, fl) — ^ 

(Fig. 1^). Here the counter part of the entropy-intruder is a region of multiple valued e(/3), this is the case for /3 
between 15 and 20. The latent heat at fixed f2, qe{^) decreases for < < 12 and is a constant for O > 12. There is 




FIG. 3. Inverse temperature /3(e, fl) surface. The mesh lines are at constant e or constant il. The intruder in S at fixed Q 
corresponds here to a multiple energy value for a given /3 and Jl, e.g. /3(e, f2 = 0) = 20 has three solutions ei ~ 0, £2 ~ —1 and 
£3 ~ —6. The thick line is close to the projection of the /3 = 00 isotherm; (3 is not defined in the forbidden region. 
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no critical value of fi, above which S{e) is concave for all e, i.e. there is a first order phase transition of all values 
of ri. In another model for a self-gravitating system such a fie was reported ||2^, but not in the one presented in [Q. 
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FIG. 4. Microcanonical angular velocity to (e, Q) surface. The mesh lines are at constant e or constant il. uj is not defined 
in the forbidden region. At high energy ui oc y/Tl', Near the ground states uj shows a richer non-monotonic behavior with peaks 
and troughs for small fl and has a nearly constant value for large Q (see text). 

On Fig. ^ we have plotted the microcanonical angular velocity w as a function of Q and e. As a direct consequence 
of Eq. (p|) UJ tends to zero with fi, and at high energy uj is proportional to (oc L). For low energies and Q < 12, 
UJ exhibits some structures with peaks and troughs, in another words at fixed e, uj is not necessarily an increasing 
bmction of fl. At high (O > 12) and near the ground states uj is almost a c onstan t. All these structures can be 
understood in terms of mass distributions which influence uj through / (see Sec. [II B). 



B. Mass distribution 



In order to understand the origin of the structures seen in the different microcanonical quantities (S', /3, w, . . .) we 
have to have a closer look at the spatial configurations, i.e. at the mass distributions. One observable we have studied 
is the mass density p (see Eq. and As the Hamiltonian of the system is rotationally invariant, the mean value 
of p can only be a function of r, the distance from the center of coordinates, although other observables might show 
a breaking of the rotational symmetry (see below) . 

On Fig. H p is plotted for different energies and for O = and = 4. For = (Fig. ||.a) we recover the classical 
case (when E is the only fixed parameter); at high energy the system is in a homogeneous gas phase (flat p), when 
the energy decreases the system undergoes a phase transition and eventually ends up in a collapse phase where a 
majority of particles are in a cluster near the center of coordinates (p peaked at r = 0). For fl ^ (Fig. |^.b) the 
situation is (a little bit) different. At high energy we recover the homogeneous gas phase. But at low energy the 
system cannot collapse entirely at the center of mass. This is due to the rotational energy e^ot = j in Eq. (|^); if the 
system contracts at the center the inertial momentum / will tend to zero and therefore erot will diverge leading to a 
negative remaining energy e^. So depending on the value of Q the main cluster will eject a certain amount of particles 
in order to increase I. Near the ground state these "free" particles will eventually collapse to form a second cluster 
in order to decrease the potential energy (j>. Due to the conservation of the center of mass, the position of the biggest 
cluster will be shifted from the center by a certain amount (see Fig. ^.b at e = —5). At low Q one particle will be 
ejected, with increasing fi the number of ejected particles raises and this process stops when two equal-size clusters 
are formed. This explains the discreteness of the peaks in uj (Fig. Q), the increase of the ground state energy eg{n); 
because the potential energy of a single cluster of size N is smaller than the one of two well separated clusters. At 
high Q > 12 the system undergoes a phase transition from a gas phase to a collapse phase with two equal size clusters 
close to the boundary. From one value of = fii > 12 to another one > the whole entropy curve at fixed 
angular momentum is simply shifted along the energy axis, i.e. S (e, ili) sa 5 (e + , i}2)- So the ground state 

energy eg{Q) at high f2 is almost on a line of equation eg + ^ + cjjg ^ 0, where and cjjg are the rotational energy 
and the potential energy of 2 clusters of size N/2 at radius r = 1 respectively. This monotonic behavior has already 
be mentioned for all the thermodynamical variables S (Fig. ||), (3 (Fig. ||), uj (Fig. ||). . . 
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FIG. 5. Density as a function of the distance from the center r at different values of energy e and angular velocity fl (arbitrary 
units). At high energy and for all fl the density is flat; the system is in the homogeneous gas phase. Near the ground state the 
density shows one peak for f2 = (a) and two peaks for £7 > (b), which correspond respectively to a one cluster and to a two 
clusters phase surrounded by some gas (see text). 



As already mentioned p is only a function of r and it can not be used to infer the angular distribution of the 
particles, i.e. there is not enough information to say if a peak in p at ro ^ corresponds to one or many clusters or to 
a uniform distribution of the particles (ring) lying on a circle of radius ro. However at least at very low energy a many 
clusters (more than two) configuration is very unlikely and will not contribute to the average values for reasons linked 
to the weight W{r) = -^Sr' For simplicity let us assume that there is only one strong peak in p at r = tq ^ 0. 
Since the center of mass is fixed this can not be the signature of a 1-cluster system. At least 2-clusters lying on 
a circle of radius ro are needed. All these n-clusters systems have the same rotational energy j = "/^j but their 
corresponding potential energy 0„ will differ. For example, with a = 0.05, j'q ~ 0.5 and N = 24: the ratio of potential 
energy is ^ ~ 1.7. So at low energy, the remaining energy e,- corresponding to a 2-clusters system will be much 
larger than the 3-clusters' one, leading to a huge difference in the weight W(r). So at low energy and for SI ^ the 
2-clusters case is dominant. At higher energies, the term Bg{I, 4>) in Eq. ( pi|) can compensate the difference in the 
weight W{r) and allow many clusters configurations and eventually at high energy a complete random configuration 
on the ring of radius ro will dominate the average mass distribution. 

We can check this argument by studying another observable, for example the normalized distance distribution P{d), 
i.e. the density of probability that the distance between two given particles is d. To probe the information given 
by P{d) we have estimated it for four simple mass distributions: (a) 2-clusters, (b) 3-clusters, (c) ring, (d) uniform 
random distribution. For (a), (b) and (c) the particle were put on a circle of radius ro — 0.5, and then randomly 
shifted several times (in order to give a spatial extension to these idealized initial configurations) ; finally the ^'^'^^^'^ 
distances are recorded for all realizations and averaged. Figure || shows the average of P{d) over 1000 realizations. 
Note that the density distribution p(r) is by construction exactly the same for the three first cases, i.e. strongly 
peaked at ro with a width of about 0.5. The latter value depends on the shift one applies on the initial idealized 
spatial configurations. 

As one can see on Fig. although the density distribution is the same for (a), (b) and (c), P{d) gives some new 
insight on the mass distribution: 

(a) There are two peaks, one at small d which corresponds to a clusterisation and another one at r = 1 = 2ro; this 
is exactly the distance between the two clusters (more precisely between their center of mass). The area under 
the small d peak is similar to the one under the large d peak, indeed the number of short distance pairs and 
the number of pairs with d ~ 1 are both about ^ . Moreover the widths of the first and second peaks are (as 
expected) ~ 0.5 and ~ 1 = 2 * 0.5 respectively. 

(b) There are again two peaks one at small d and another at d ~ 0.8 < 1 and their respective widths are similar to 
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FIG. 6. Average of P{d) the distance distribution for different simulated spatial configurations. See text. 



the ones in (a) . The large d is compatible with the length of one side of the equilateral triangle on top of which 
the 3-clusters mass distribution has been built. This time the area under the large d peak is larger than the 
one under the short d peak, since the number of short distance pairs is about whereas the number of pairs 
with d~ 0.8 is 

(c) In the ring a trace of the two peaks still exists but they are not well separated because a lot of intermediate 
distances are compatible with this model. 

(d) When the particles arc uniformly distributed Pid) is a binomial-like distribution. 



£i=0 


E=100 




£2=0 


E=0 




- £2=0 


E=-5 




-- £2=4 


E=-5 






FIG. 7. Distance distribution P(d) for the gravitational case for different values of (e, fi). At high energy P{d) corresponds 
to a random distribution (see Fig. 0). For = and at low energy, P{d) has one peak at d ~ 0; almost all particles are very 
close from each other, and there is a single cluster collapse phase. For Q there are two peaks at low energy: one at very 
small d which is a sign of clusterisation and another peak at large d which signals multiple clusterisations; in fact there are two 
clusters (see text). 

We have also estimated P{d) for our gravitational system, as shown on Fig 0. At high energy, P{d) corresponds 
to the randomly distributed case (see Fig. At low energy with = 0, P{d) has only one peak at d = 0, this 
corresponds clearly to a single cluster case surrounded by some gas. For ^ and at low energy (in Fig. Q e = — 5 
and = 4), there are two well separated peaks, one at small do = and the other at c?i ~ 1.1. The peaks imply the 
presence of at least two clusters, however the fact that the widths of the peaks are small excludes a large number of 
clusters and even more the ring case (see Fig. ^). Now we can combine these informations with the ones obtained by 
studying p(r) (see Fig. ||). For e = — 5 and fl = A, p has two peaks at ri ~ 0.15 and r2 ~ 1. All in all, this means 
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that there are, in the mean, two clusters rotating around the center of mass. The distance between these clusters is 
ri+r2 — 1.15 ~ di. Their mass ratio is — = — ~ 0.15. Since we know the total mass mi +m2 = 20 we set mi ~ 17 
and m2 ~ 3. 

The distance distribution can be of great help to identify the mass distributions at low energies. However at the 
transition regions since there is a superposition of different types of mass distributions the knowledge p and P{d) 
is not sufficient and therefore of no help if we want to study for example the "fractality" of the mass distribution 
as it has been discussed in other self-gravitating systems p0|-p2t (not mentioning the fact that the actual number 
of particles in the presented numerical applications is too small) , and further work is needed to get a more detailed 
picture. 

At very low energy, near the ground state at least one of the clusters (the smallest) is very close to the boundary. 
There the assumption of a small evaporation rate made in Sec. | does not hold. 



C. Phase diagram 



a 10 




Xf>0. X2<0 
^1 Xf>0. X2>0 
\/'/'\ G region 



FIG. 8. sign (Ai(e, Q.)) as defined in Sec. [IE . The white regions correspond to Ai < and therefore to a positive determinant 
Ds = A1A2 > of the entropy curvature, these are pure phase regions. The gray region corresponds to Ai > and Ds < 
and the dark gray ones also to Ai > but Ds > 0; Ai > defines first order phase transition regions (see text). Points in 
G correspond to local maxima (minima) in Eq. ( p^ ) if Ai < (Ai > 0) but not to global maxima (minima); points outside G 
correspond to global maxima in Eq. There is a one-to-one mapping between the microcanonical ensemble and the GBE 

only in the white region (see text). S is not defined in the forbidden region, here in light gray. Note that (a) the points at 
57 = and low energies e < —7 are not included in G, (b) the high energy limit of G is only known approximatively. 



On Fig. g we have plotted sign (Ai(e, fi)) as defined in Sec. II B. This plot can be taken as the phase diagram of 
the self-gravitating system at fixed e and fl. The white regions correspond to pure phases (Ai < ). At high energy 
there is a homogeneous gas phase and at low energy there are several pure collapse phases with one (fi = 0) or tw o 
(O ^ 0) clusters. The diff'erent 2-clusters phases are characterized by the relative size of their clusters (see Sec. Ill B| ) . 
These regions are separated by first order phase transition regions where Ai > (gray in Fig. There are even two 
regions (dark gray) where the entropy S" is a convex function of e and uj; i.e. all the eigenvalues of Js are positive 
(Ai > and Ds > 0). These two regions are quite stable with respect to the number of particles (at least for N < 50). 
Their surface increase slightly with the number of particles N. The orientation of Vi the eigenvector associated with 
Ai (defined as the largest eigenvalue of Js) is not yet known in detail for all (e, fi); however we can already state that 
at "high" energy Vi is almost parallel to the energy axes (phase transition in the e direction) and should be parallel 
to the ground state at very low energy. The overall structure of the collapse phases matches the one of the angular 
velocity to (see Fig. |4|): roughly, the peaks in u correspond to pure phases while the valleys between these peaks belong 
to the first order phase transition region. Due to a lack of precision we are yet unable to ensure that the two isolated 
pure phases at e « —4, fi « 1 and O « 7 are really isolated and not linked to a pure phase at lower energies. 

As already mentioned, unlike in the model presented by Laliena in pst there is no critical angular momentum Lc 
above which the first order phase transition vanishes giving rise to a second order phase transition at Lc- Nevertheless 
this does not exclude second order phase transit ion (critical point) at all. They are defined in the microcanonical 
ensemble by: (i) Ai — 0; (ii) VAi • vi = (see Sec IIB). On Fig. || (just like on Fig. ^) regions where Ai < (> 0) are 
in white (gray). The condition (i) is simply achieved at the boarder between the gray and the white regions. The thick 
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FIG. 9. Locus of second order phase transitions (see text). 



lines on Fig. |^ correspond to condition (ii) . Second order phase transition are located at the crossing points points of 
the thick lines and the boarders. One immediately sees that there are several critical points. However there are not all 
of (astro-)physical interest since most of them are close to the ground states line or at very high angular momentum 
where the small evaporation rate assumption is not valid. Nevertheless there are two points one at (e, 51) « (—0.5, 1) 
and another one at (e, f2) « (—0, 4) where this assumption is valid and therefore their deserve further investigations 
and especially regarding their corresponding mass distributions. 



D. Loss of information in the {(3,13^) ensemble (GBE) 

In a recent paper Klinko and Miller have studied another model for rotating self-gravitating systems jH]. They 
introduced the canonical analogous of the (e, fi) ensemble namely the (/3, 7/3) ensemble (GBE), see Eq. (^) and (||). 
Once one defines the generalized microcanonical ensemble (ME), it is straightforward to introduce and study the 
system in its conjugate ensemble, the generalized canonical ensemble (CE). Note that "generalized" means in CE "as 
a function of all the intensive variables" and, in ME "as a function of all extensive parameters" . There are mainly 
two reasons invoked to study a system in CE instead of in ME 

1. Performing the computations in CE arc in most cases much easier than in ME. Here CE can be seen as a trick . 
However there is a priori no reason for the results to be equivalent. Indeed the strict equivalence of the ensemble 
is only achieved at the thermodynamical limit except in the first order phase transition regions p^j3^ , p^ . In a 
weaker sense, needed for Small systems, ME and CE can only be equivalent when Xi{X) < and far from any 
first order phase transition regions (see below). Moreover recent progress in computer performance enables one 
to perform now numerical experiments within the ME for increasingly complex systems. 

2. The studied system is not isolated but is in contact with a "heat bath" and exchange amounts of X with it. 
Then, obviously ME does not provide a suitable description, and CE might be eligible. Here CE describes a 
different physical system than the one in ME. At the thermodynamical limit (if it exists) the ensembles are again 
equivalent except at first order phase transitions. However, the CE description is valid only if the Hamiltonian 
of interaction Timt is small comparing to the ones of the system Tisys and the heat bath Tiht- This condition is 
usually fulfilled at the thermodynamical limit (if it exists), but for Small ones Timt can hardly be small at least 
comparing to Hsys and this can lead to dramatic effects [ |5^ , |55[ |. In this case a better description would be the 
ME of the system, its "heat bath" and their interactions. 

For now on we consider in this paper the cases when CE is used as a mathematical trick, and we focus and the 
amount of information lost from ME to CE. 
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The link between ME and CE is given by Laplace transform (using the notations introduced in Sec. [IB) 



/ dXe-^-^+^W, (25) 



where Z is the partition sum of the CE; x — {z^, . . . are the intensive variables associated with X and defined 

by 

^ - OX^^ (26) 

for i = 1, ... ,7V. The ME at Xq is equivalent to the CE at xq = x{Xo) if the integrand in Eq. (||) e^^-^'^^o^ and 
therefore /(X, a;o) = —X ■ xq + S{X) has a global maximum at Xq. This condition is violated when Xi{Xo) > (this 
is the basic idea behind the definition of phase transitions in Small systems, see p^). So in practice all informations 
contained in points where Ai > are lost after the Laplace transform (^^ in CE. 

Now we are left with the points Xq characterized by Xi{Xo) < 0; this relation implies only that f{X,xo) has a 
local maximum at Xq but not that it is a global one; Ai < is a necessary but not sufficient condition. In Fig. |l^ 
we have illustrated this point with a trivial one-dimensional example; all the points below the Maxwell line do not 
correspond to global maxima and their information content is smeared out and, in practice, lost in CE. 




FIG. 10. Schematic entropy curve S as a function of one extensive parameter X (solid line) and the Maxwell line (dashed 
line). A and D are the end points of the Maxwell construction. At B and C the largest curvature vanishes, i.e. Ai = 0. All the 
points below the Maxwell line do not not correspond to a global maximum in the Laplace transform Eq. (^). This includes 



points where Ai 



> 0, but also points where the entropy is concave, for X G [Xa,Xb] and X £ [Xa,Xb\, Ai(X) < 0. 



From the last remark we see that a way to check if a point corresponds to a global maximum in Eq. ( p5| ) is to study 
what would be a generalization of the Maxwell line in A/" > 1 dimensions. Work is in progress in this direction and 
the results will be presented elsewhere, but we can already state that this task is, to some extent, similar to the one 
of building the convex hull of a set of points in M dimensions, or more specifically the convex hull of the entropy S. 

However, if one needs qualitative results, in two dimensions the task can be rather easily solved in another way. 
Let us define p{X,Xq) the tangent plane to S[X) at Xq] its equation is p{X,Xq) = —X ■ xq + S{Xq) + Xq ■ xq. If 
/(AT, Xq) has a global maximum at Xq then 



p{X,Xq) > Six) 



(27) 



for all X ^ Xq and p{Xq, Xq) = S{Xq) by definition. 

In the case of the gravitational model presented here X = {e, fi} and x = {/?, 7/3}. Now if one inspects the entropy 
surface S{X) (see Fig. it is clear that condition (|27|) is not satisfied for all the points in the regio n filled with 



dashed lines (G) in Fig. g^. This is due to the concavity of the energy ground state eg{n) (see Sec. Ill A). G includes 



all the two-clusters collapsed phases, the first and second order phase transitions (except for lo — Q — j = 0). All 
the information contained in G is smeared out through the Laplace transform (ESI) in GBE and, in practice, lost. 
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The fact that GBE misses all the two-clusters collapse phases would be already enough to disqualify it as being 
a good approximation (mathematical trick) of the ME. But, furthermore, if one studies more carefully /(X, xq), 
Xq e G; one will notice that (a) there is one local maximum at = and (b) there is no maximum for high fi: in 
the direction of increasing Q at low energy, f(X,xa) is a never ending increasing function, i.e. /(X^xq) has no global 
maximum for X G G (see Fig. |ll|). Therefore the integral in Eq. ( p5| ) diverges for all xq, Xq G G. In other words the 
GBE, in our model, is not defined for high (3 and 7 7^ {uj ^ Q). 




FIG. 11. /(e, n, /3o, 70) = —e/3o— f2/3o7o+5'(e, if as a function of e and fi, where if is an arbitrary constant; = /3(eo,no); 
f^o = f2(eo,no); Xq = (eo,r2o) = (—3,5) ^ (f2o,7o) ~ (27.9,-0.196). The mesh lines are at constant e or constant Q,. As 
expected / has a saddle point at Xq « (/9o, 7o) since m(Xo) > and Ds{Xo) < (see text and Fig. / has a global maximum 
at n = and e « —7, but one sees that it is an monotonically increasing function for increasing Q, and e{Q.) ~ ^g{i^) + 2. 
Therefore the integral in Eq. (p3) diverges and the (/3, 7/?) ensemble is not defined for (/3o,7o). 



This example shows how dramatic can be the information loss if one studies an isolated system in the CE. 



IV. SUMMARY AND DISCUSSION 



A. Summary 

The aim of this paper is to present the results of a study of a self-gravitating system of N classical particles on a 
disk of radius R for which two extensive variables, namely total energy E and total angular momentum L are fixed. 
The microcanonical entropy S can be written as the logarithm of an integral over the spatial configurations (^ and 
(^ . The conservation of L implies that the mean value of the linear momentum p of a particle at a radial distance 



r is proportional to L and r (13); in the mean the system rotates like a solid body, i.e. the mean angular velocity of 
a particle does not depend on its position (p7|). The dispersion of p is broader than the one it had if L would not 
be conserved and depends on the radial position (^8|). In order to integrate e"^ we write it as the folding product of 
a "background" function Bg and the weight associated with the remaining energy (^). Once a numerical estimate 
of Bg is obtained, S and its derivatives, such as the inverse temperature /?, can be computed. The entropy surface 
shows an intruder in the energy direction which signals a first order phase transition with negative specific heat (Fig. ^ 
and Fig. ||). Contrary to another model there is no critical value of L above which this transition is no longer 
present, nevertheless there are several points of the parametric space (E, L) . At high energy the angular velocity lo is 
a simple increasing function of L (Eq. (^) and Fig. ^ , but at low energy, near the ground states the relation between 
L and to becomes non-trivial. All these peculiarities can be understood if we study the mass distribution. We use 
two observables: the density of mass p as a function of the radial distance; and P{d) the density probability that the 
distance between two particles is d. At fixed L and at high energy the particles are randomly distributed over the 
disk (gas phase); if the energy decreases, at some threshold the system undergoes a phase transition of first order to 
a collapse phase. If L = there is a single cluster at the center of mass, if L ^ then there are two clusters rotating 
around the center of mass. Their relative mass ^ with mi > TO2 is very large for low L, it decreases as L raises and 
eventually becomes 1 for L larger than a certain threshold. Clearly these two clusters phases were not reported in 
previous work because of the usual a priori assumption of rotaional symmetry. 
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Phase transitions (and phases) can be defined unambiguously in the microcanonical ensemble as a function of the 
local topology of S. If all eigenvalues of Js, the Jacobian of S are negative then the system is in a pure phase. If 
(at least) one eigenvalue is positive, say Ai, then there is a first order phase transition in the direction of vi, the 
eigenvector associated with Ai. A critical behavior occurs when one (or many) eigenvalue vanishes on at least a second 
order region (if A^ — 0, Vv^A^ — 0). Note that one can find positive, negative and zero eigenvalues at the same point 
in the parameter space. Using these criteria, we can draw a phase diagram of the gravitational system (Fig. |8|) at 
fixed E and L: at high energy there is a pure gas phase; at low energy there are several pure collapse phases, there 
is one cluster at L = and there are two clusters for L 0; these phases are separated by first order phase transition 
regions. There are also several second order phase transitions; two of them are located at relativley high energies 
{E « 0) and therefore may be of astrophysical importance. 

Studying an isolated system using the canonical ensemble (CE) can be very misleading, since there is a massive 
loss of information from the correct microcanonical ensemble (ME) description to the CE's one (if there are phase 
transitions in ME). In fact for the gravitational model, CE cannot be sensitive to the two clusters phases and all t he 



phase transitions and it is not defined for some values of the intensive parameters that exist in ME (see Sec. [II D) 



B. Discussion 

Of course we have just presented an equilihirum statistical model that may help to understand the physics of 
globular clusters or collapsing molecular clouds and the results should be interpreted with caution especial in the case 
of star formation. A lot of "ingredients" are missing in order to have a complete picture of the formation of multiple 
stars systems and planetary systems, for instance the magnetic field p6| , ^ , or the presence of vortices js^]. 

Phases and phase transitions can be well defined in the microcanonical ensemble (ME) without invoking the ther- 
modynamical limit by probing the curvature of the entropy surface. There still exists open problems. One of them is 



the scaling of a first order transition: in Fig. 10 in a canonical sense the phase transition occurs from Xa to Xo (if 
X is the energy E then Xo — Xa is the transition latent heat), but for if one uses the ME definition the transition 
occurs only from Xb to Xc- At the thermodynamical limit, if it exists, such discrepancy should disappear. In another 
context (a model of a first order liquid gas transition of finite-size Sodium clusters) we could show that 1 
when the system size goes to infinity. 

The ME definition of phase transitions offers a richer view of physical systems and phase transitions. Again on 
Fig. CE is not sensitive to everything that could happen under the Maxwell line: there is always one transition. 
On the contrary, in ME there could be many phase transitions between Xa and Xd- For example if there is a small 
positive curvature bump between Xb and Xc there would be two transitions in the microcanonical sense but still 
one in the canonical sense. 
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APPENDIX A: 

Let us compute {Pk)q^ the average momentum of particle k at fixed position (for simplicity we will set fc = 1). The 
a-component of {Pi)q-^ is 



(Pi] 



I 




pfs {E - n) 6^ (E. p.) s (E. qiXPi- l) 5^ (E. q^) 




)S{E-H)S^ (E. P^) S (E. Q^ X p. - L) 5^ (E, Q.) 



^ ^ (AI) 



I{Ut2dq.)W{r)5^ (E.g.) 

where = J (Ylidpi) PiS {E — H) {J2iPi) ^ {J2iQi ^Pi~L), and W{r) is the microcanonical weight at fixed 
spatial configuration r, its value is W{E,L,r) = C^E^~^^^ (see (^)). The outline of the derivation of is the 
same as in [psi for W, and we get after some algebra 



15 



JV-5/2 



(5=1 



where e is the antisymmetric tensor of rank 2. Using ( A2) in (Al) we get finaUy 



/ (nf=2 rfg.) ^-^/^ ELi rfesc. Er"'5^ (E, 9.) 
= Lmi{r^)r^^r{esa. 



(A2) 



(A3) 



FinaUy 



{Pi)ri = Lmi{I ^)rk^rfesaea, 

6, a 



where is the a-component unit vector. 

(Pk)Qk '^^^ derived in a similar way, and we get 



2mi 



{N-5/2){Er-^)q^ 

mkh 

IiN-5/2){Er-^)q^ 



M ) 



(A4) 
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